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Evaluation of a vortex-based subgrid 
stress model using DNS databases 

By Ashish Misra 1 AND Thomas S. Lund 2 


The performance of a subgrid stress (SGS) model for large-eddy simulation (LES) 
developed by Misra &z Pullin (1996) is studied for forced and decaying isotropic 
turbulence on a 32 3 grid. The physical viability of the model assumptions are 
tested using DNS databases. The results from LES of forced turbulence at Tay- 
lor Reynolds number R\ ~ 90 are compared with filtered DNS fields. Probability 
density functions (pdfs) of the subgrid energy transfer, total dissipation, and the 
stretch of the subgrid vorticity by the resolved velocity- gradient tensor show rea- 
sonable agreement with the DNS data. The model is also tested in LES of decaying 
isotropic turbulence where it correctly predicts the decay rate and energy spectra 
measured by Comte- Bellot &; Corrsin (1971). 


1. Introduction 

The aim of this study is to use DNS data of isotropic turbulence to evaluate the 
performance of a new subgrid- stress model. LES is performed for both forced and 
decaying turbulence, and statistics are compared with appropriately filtered DNS 
fields and with some experimental results. The SGS stresses are calculated from 
a structural model of the subgrid vorticity proposed by Pullin &: Saffman (1994), 
henceforth PS. This model has some similarity to the eddy-axis structure model for 
one-point closure proposed independently by Reynolds &; Kassinos (1995). PS as- 
sume that the subgrid structure consists of an ensemble of straight stretched vortex 
structures each with an arbitrary internal vorticity distribution. Some support for 
this type of structure of the fine scales is provided by the observed tendency, in sev- 
eral numerical simulations, for the alignment between the vorticity vector and the 
eigenvector corresponding to the algebraically intermediate value of the principal 
rate-of-strain. Moreover, stretched- vortices have been used to make quantitative 
predictions for a range of fine-scale turbulence properties (Lundgren 1982, Pullin 
&: Saffman, 1993, 1994). Misra Sz Pullin (1996) have examined and implemented 
several different versions of the locally anisotropic model of PS. In this report we 
will discuss one such model and examine its performance when measured against 
experiment and filtered DNS. 

1 Graduate Aeronautical Laboratory, California Institute of Technology, Pasadena CA 91125 

2 Center for Turbulence Research 


360 


Ashish Misra & Thomas S. Lund 


2. Vortex orientation model 

In the PS model, the orientations of the structures are given by a pdf of the distri- 
bution of Euler- angles describing the transformation from laboratory to structure- 
fixed axes. The Reynolds stresses are proportional to the turbulent kinetic energy 
of the vortex collection times a tensor-moment of the pdf. For example, a delta- 
function pdf in which all vortices in a subgrid domain have a common direction 
described by the unit vector e gives Reynolds stresses, 
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where k c is the cutoff wavenumber and the subgrid kinetic energy is given by K = 
E(k)dk. PS proposed a simplified version of the model using a rapid-distortion- 
like approximation in which the orientation of the subgrid vortices have a two- 
delta function pdf defined by the eigenvectors corresponding to the largest two 
eigenvalues of the resolved rate-of-strain tensor S t) . It can be shown that such 
a model can not produce backscatter. In this report we study the performance 
of a slight variant of their original model wherein there is alignment between the 
eigenvector corresponding to the maximum eigenvalue of the rate-of-strain tensor, 
e^, and the resolved vorticity vector, Co. The Reynolds stresses are then given by 
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where // is the fraction of structures aligned with the maximally extensive eigen- 
vector and is the unit vector along Co. As partial justification for (2) we remark 
that one should expect complete alignment with Co in the DNS limit. We currently 
take fj, = 0.5. In order to calculate the subgrid energy, A , a local balance be- 
tween production by the resolved scales, and the sum of subgrid and resolved scale 
dissipation is assumed. When coupled to an assumed Kolmogorov subgrid energy 
spectrum produced by the (unknown) internal structure of the vortices with a cut- 
off at ki] = 1, where 7] = (i/ 3 /^) 1 ^ 4 is the local Kolmogorov length, an equation 
sufficient to determine the dissipation is obtained. This is given by 
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where A'o is the Kolmogorov constant. When the model parameter A () is specified, 
(3) can be solved for the total dissipation e and the subgrid energy determined from 
the Kolmogorov spectrum. This gives closure. 

Equation (3) has the dimensions of L 2 T~ 3 ; we therefore divide (3) by k*v 3 which 
results in two non-dimensional parameters (see Misra A.' Pullin for details), 
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where S\ represents the resolved scale dissipation while £2 represents the stretch 
experienced by the subgrid vortices by the resolved velocity-gradient tensor. The 
SGS dissipation may be written as, 


= - e 3 ,e 3j ) + \(S tJ - A', 
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where S 3tr is the component of Sij aligned with the vortex. Hence backscatter, 
defined by € sg9 < 0 occurs whenever S 2 > 0 - the subgrid vortices are being com- 
pressed on the average - while 52 < 0 - the vortices are axially stretched - gives 
cascade. 

This model has been implemented for both forced and decaying box turbulence 
by Misra & Pullin. They examine several alternative scenarios for determining the 
instantaneous orientations of vortices in a given cell. 

3. Results and discussions 

The incompressible, filtered Navier-Stokes equations are solved in a 32 3 box, with 
and without forcing, using periodic boundary conditions in all three directions. A 
Fourier-Galerkin pseudo-spectral method is used with c 3/2 dealiasing rule 1 for the 
non-linear terms, i.e. 32 Fourier modes in each direction were advanced in time, the 
computation of the non-linear terms were done using 48 modes in each direction. 
A second order explicit Runge-Kutta scheme is used for time advancement. 

3.1 Decaying turbulence 

We study decaying isotropic turbulence in order to compare our results to the 
experiment of Comte-Bellot and Corrsin. They measured the energy spectrum at 
three downstream locations in grid turbulence. One can relate this to decaying 
isotropic turbulence by invoking the Taylor approximation. We mimic their exper- 
iment by studying turbulence in a cubical box with periodic boundary conditions. 
In a frame of reference moving with the mean flow speed, 

Jo V(j') 

where x is the downstream distance from the grid and U(x) is the mean flow ve- 
locity over the cross-section of the tunnel. We have non-dimensionalized the ex- 
perimental data by the following characteristic velocity, length and time scales: 

U re f = yj 3Aq 2 /2, L re j = L/2 tt and t re f = L re j jUrtf * In their experiments the 

velocity fluctuation at the first measuring station is \/Uq 2 = 22.2 cm/s, the free- 
stream speed is Uqq = 10 3 cm/s and the spacing of the turbulence generating mesh 
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FIGURE 1. Decay of resolved ( ) and subgrid ( ) kinetic energy. No 

model: . Experiment: Resolved, ■; subgrid, ▲ . 



Figure 2. Time evolution of spectra in decaying turbulence. Experiment: □ , 
t =42; A , t =98; o , t— 171. Simulation at corresponding times: . 

is M = 5.08 cm. The size of the computational box, L — 11 M, was chosen to 
contain roughly four integral scales. The times at the three stations were measured 
in terms of Uoot/M. The initial Taylor Reynolds number is R\ ~ 80. In order 
to compare the resolved and the subgrid part of the turbulent energy produced by 
the computation, the measured spectra have been integrated over the relevant scale 
ranges. 

Figure 1 shows the decay of the resolved energy with time. The LES gives good 
agreement with experiment. The dotted line is the result of running the simulation 
with the model switched off. It is evident that the model plays an important role in 
predicting the correct decay of the kinetic energy. Aside from the parameter // (set 
to 0.5), the model requires a value of the Kolmogorov constant. While acceptable 
results were obtained with the standard value of Jv 0 = 1.5, a slight improvement 
was observed when higher values were used. The results in Fig. 1 were obtained 
with A' 0 = 1.8. While this value is on the edge of the uncertainty band from 
experimental measurements (Sreenivasan 1995), it is well within the predicted range 
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from numerical simulation; Jimenez et al report a value of A'o « 2 in a 256 3 DNS 
of isotropic turbulence at R\ — 94.1. It is possible that elevated values of the 
Kolmogorov constant in numerical simulations is a low Reynolds number effect. If 
this is the case, it stands to reason that a larger value should be used in the present 
simulations which are also performed at relatively low Reynolds number. 

The decay of the subgrid energy with time is also shown in Fig. 1. Note that 
the subgrid energy is obtained from the model without the solution of additional 
field equations. The subgrid energy is a quantity derived from a knowledge of 
the resolved field and the chosen subgrid energy spectrum; it therefore can not be 
initialized to match the experimental value. Figure 2 shows a plot of the resolved 
energy spectra with the measurements at the initial time and then at the two later 
instants. The initial spectrum is generated to match the experimental data, while 
the later two curves are the predicted spectra arising from the LES calculation. 
Figure 1 gives the area under the curve of Fig. 2 at the three time instants, over 
the resolved range of scales. 


8.2 Forced turbulence 


Forcing is achieved by exciting low wavenumbers such that the total energy in- 
jection rate is constant in time. A certain selected number of Fourier modes are 
chosen from a wavenumber shell |k| = k 0 . The Fourier coefficient of the forcing 
term is then written as, 




6 U* 

N ~ 2 

|U*| 
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for all modes in the specified shell. The above choice of ensures that the energy 

injection rate, ^ f* • U*, is a constant and equal to 6. We have chosen Jfc 0 —2 ,N = 
20, and 8 = 0.1 for all the runs. (See Carati et al. 1995) 

Simulations with forcing were performed such that a statistical steady state is 
reached when statistics are collected. Results in this report are presented for Taylor 
Reynolds number, R\ ~ 85 in order to make comparisons with 128 3 DNS results at 
approximately the same Taylor Reynolds number. 

Figure 3a-b shows scatter plots of Si versus S 2 from the LES as well as from 
filtered DNS data. These plots show the intensity of the vortex stretch as a function 
of the resolved dissipation rate. Notice that the DNS data displays a significant 
fraction of points with positive stretch parameter, S 2 , (backscatter) whereas the 
model rarely predicts these events. Quantitatively, the DNS shows roughly 30% 
backscatter, which is consistent with previous measurements (Piomelli et al. 1991). 
In contrast, our model yields only ~ 3% backscatter. While there is clearly a large 
discrepancy in the prediction of backscatter in Fig. 3, it should be noted that the 
percentage of backscatter can be controlled through the parameter ji in Eq. 2. \i = 1 
corresponds to complete alignment with 63 and results in no backscatter, while 
^ = 0 corresponds to complete alignment with cl? and gives about 40% backscatter. 
When 0 < /i < 0.4, the decay of the kinetic energy appears correct, but the decay of 
the energy spectra is somewhat unsatisfactory with a trend towards flattening of the 



364 


Ashish Misra & Thomas S . Lund 



Figure 3. Scatter plot of Si and S 2 in filtered DNS (a) and LES (b). 


spectrum at later times. In the range 0.4 < fi < 1, there is a general insensitivity 
to leading to correct statistics and diminishing backscatter. The performance of 
the model for all values of fi in forced turbulence is satisfactory. Presently we show 
results for fx = 0.5 which are typical of the behavior of the decay of the resolved 
energy spectrum and of the resolved and subgrid energy for 0.4 < fi < 1.0. 

An interesting feature of Fig. 3 is that all points lie within a bounding parabola. 
An estimate based on a locally two- dimensional ‘maximum stretch 1 scenario for S tJ 
gives a bounding parabola Si = 16 S% . We find however that Si = 12 S% gives a 
slightly better boundary and so this curve is displayed in the figure. The importance 
of backscatter has been a question of debate though there is some evidence in the 
literature of its importance in wall-bounded flows. The backscatter property of the 
model is also illustrated in Fig. 4a, which shows a pdf of the ‘stretch 1 , that part 
of the velocity- gradient tensor which stretches the subgrid vorticity. The stretch is 
suitably normalized by \/e/u. While the two curves peak at approximately the same 
location, the LES shows predominant stretching. Figure 4b is a plot of the pdf of the 
subgrid energy transfer, e sgs /t] points on the left of the origin exhibit backscatter. 
It is clear from the figure that the LES does well in the cascade region. The pdf 
of the total dissipation logio (e/e) is displayed in Fig. 4c. The total dissipation e is 
a positive definite quantity by construction as defined in (3). The distribution of e 
appears to be approximately log-normal. 


4. Concluding remarks 

The behavior of an SGS model for LES has been tested against filtered DNS fields 
at similar Reynolds numbers. The model is stable and appears to produce a good 
quantitative description of the resolved flow and of the subgrid energy. It shows 
the right decay rate and gives good agreement with the experiments of Comte- 
Bellot and Corrsin. The model seems to work well for forced turbulence. Misra 
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FIGURE 4. Pdf of: (a) stretch, (b) subgrid energy transfer, (c) total dissipation. 
LES, DNS. 

& Pullin demonstrate that it produces negligible SGS dissipation in the limit of 
fully resolved flow. A small amount of backscatter is produced by the version of 
the model tested presently albeit not as much as is indicated by the filtered DNS 
field. An adjustable parameter in the model is the Kolmogorov constant. The 
value presently used is within the bounds of experimental values. The model is 
about 25% more expensive in CPU time than the simple Smagorinsky model. Some 
subgrid model features show qualitative but not strong quantitative agreement with 
the equivalent quantities from filtered DNS fields at a similar Reynolds number. It 
remains to be seen how well the model will function for free-shear or wall-bounded 
flows. Future work will aim at constructing models with alternative representations 
of the subgrid vortex orientations. 
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